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Abstract. Entropy functionals arc computed for non-stationary distributions of par- 
ticles of Lorentz gas and hard disks. The distributions consisting of beams of particles are 
found to have the largest amount of entropy and entropy increase. The computations show 
exponentially monotonic increase during initial time of rapid approach to equilibrium. The 
rate of entropy increase is bounded by sums of positive Lyapounov exponents. 

1 Introduction 

The H-theorem for dynamical systems describes the approach to equi- 
librium, the irreversibility and entropy increase for deterministic evolutions. 
Suppose that a dynamical transformation T on a phase space X has some 
"equilibrium" measure //, invariant under T, i.e. n(T~ x E) = fi(E) for 
all measurable subsets E of X. Suppose also that there is some mixing 
type mechanism of the approach to equilibrium for T, i.e. there is a suf- 
ficiently large family of non-equilibrium measures v such that Ut(E) =: 
v(T~ l E) —>-t^oo n{E) for all E. Then, the H-theorem means the existence 
of a negative entropy functional S(u t ) which increases monotonically with t 
to zero, being attained only for v — \i. The existence of such functional in 
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measure-theoretical dynamical systems has been the object of several inves- 
tigations during last decades see 0- 0, [Hj [JU UH1 HE])- Here we study 
this problem for the Lorentz gas and hard disks. The dynamical and stochas- 
tic properties of the Lorentz gas in two dimensions which we consider here 
was investigated by Sinai and Bunimovich as an ergodic dynamical system 
[To! HI |3] . Other transport properties have been also studied numerically (see 
[TU]). This is a system of non interacting particles moving with constant 
velocity and being elastically reflected from periodically distributed scatter- 
ers. The scatterers are fixed disks. On account of the absence of interactions 
between particles the system is reduced to the motion of one billiard ball. 
We shall investigate the entropy increase under the effect of collisions of the 
particles with the obstacles. For this purpose, we consider the map T which 
associates to an ingoing state of a colliding particle the next ingoing collid- 
ing state. The particle moves on an infinite plane, periodically divided into 
squares of side D called " cells" , on the center of which are fixed the scatterers 
of radius a ( Fig. |TJ). The ingoing colliding state is described by an ingoing 
unitary velocity arrow at some point of the disk. To a colliding arrow Vi(Pi) 
at point Pi on the boundary of the disk the map associates the next colliding 
arrow V2(P2) according to elastic reflection law. Thus, the collision map 
does not take into account the free evolution between successive collisions. 
Let v be a non-equilibrium measure, which means that v is a non invariant 
measure approaching the equilibrium /i in the future. It is mathematically 
possible to define a non-equilibrium entropy for a family of such measures, 
using conditional expectations (i.e. a generalized averaging) relatively to the 
some remarkable partitions, namely the contracting fibers of the hyperbolic 
dynamics However, in our numerical simulations some given finite 

precision is needed, so that we consider partitions into cells with positive 
/x-measure. Here, we use slightly similar entropy functionals. Starting from 
the non-equilibrium initial distribution u, and denoting by V such partition 
formed by cells (Vi,V2, V n ) and by u t (Vi) — v o T~*(Pj), the probability 
at time t for the system to be in the cell V% and such that v{Vi) 7^ /^(Pj) for 
some i, the approach to equilibrium implies that u t (Vi) — > /i(Pi) as t — > 00 
for any i. The entropy functional will be defined by: 

s(t, u,v) = - g ut{Vi) M^j) := -n(t, v, v) (1.1) 

which we simply denote here after S(t). The H-functional (jl.lj) is maximal 
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Figure 1: The motion of the particle on a toric billiard. 



when the initial distribution is concentrated on only one cell and minimal if 
and only if v t (Vi) = fi{Vi), Vz. These properties are shown straightforwardly. 
This formula describes the relative entropy of the non-equilibrium measure 
u t with respect to fi for the observation associated to V. It coincides with the 
information theoretical concept of relative entropy of a probability vector (pi) 
with respect to another probability vector (qi) defined as follows: — In pi be- 
ing the information of the i th issue under the first distribution, — YliPi ^ n (~)) 
is equal to the average uncertainty gain of the experience (pi) relatively to 
(ft)- 

A condition under which formula shows a monotonic increase with 
respect to t is that the process u t (Vi) = v o T~ l (Vi) verifies the Chapman- 
Kolmogorov equation valid for Markov chains and other infinite memory 
chains. For a dynamical system, this condition is hardly verified for given 
partition V. However, the very rapid mixing leads to a monotonic increase of 
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the above entropy, at least during some initial stage, which can be compared 
with the relaxation stage in gas theory. 

In this paper, we will first compute the entropy increase for some remark- 
able non-equilibrium distributions over the phase space of the Sinai billiard. 
The billiard system is a hyperbolic system (with many singularity lines) and, 
in order to have a rapid mixing, we will consider initial distributions sup- 
ported by the expanding fibers. Such initial measures have been used in 
HH EE] • For the billiard the expanding fibers are well approximated by 
particles with parallel arrows velocity. We call this class of initial ensemble 
beams of particles. We first compute the entropy increase under the collision 
map for these initial distributions. We will consider finite uniform partitions 
of the phase space as explained below. The entropy functional will be de- 
fined through (jl.ljl . For this purpose, the phase space of the collision map 
is described using two angles (fl,ip), where f3 is the angle between the outer 
normal at P and the incoming arrows V(P), (3 G [0, |[, and ip G [0, n] is 
the angle between x-axis and the outer normal at P. Thus, the collision 
map induces a map: — > (Z^,^) (see Fig. H3J) and we shall first 
use a uniform partition of the (/3, ip) space. The computation shows that 
whatever is the coarsening of these partitions the entropy has the monotonic 
property in the initial stage. It is clear that, along mixing process, the initial 
distribution will spread over all cells almost reaching the equilibrium value. 
Physically, this process is directed by the strong instability, that is expressed 
by the positive Lyapounov exponent. 

We also consider the relation of the rate of increase of the entropy func- 
tionals and Lyapounov exponents of the Lorentz gas. Our computation shows 
that this relation is expressed by an inequality 

max(S(n + 1) - S(n)) = AS < ^ A; (1.2) 

Ai>0 

where the "max" is taken over n, which means that the K-S entropy is an 
upper bound of the rate of increase of this functional. 

In section 3, we shall consider another phase space and another partitions 
associated to spatial extension of the motion of the Lorentz gas. Here the 
space in which moves a particle is a large torus divided into rectangular cells, 
in the center of each cell there is one disk. Denoting the total number of cells 
by n and the number of particles initially distributed in only one region, 
by N, and following them until each executes t collisions with obstacles, we 
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compute the probability that a particle is located in the ith cell as given by: 
Number of particles in celH having made t collisions 



Pi(t) 



N 



The equi-distribution of the cells leads to take, as equilibrium measure, /i, 
-, so that this "space entropy" is defined by: 



n 



S sp (t) = -}_^p l (t)\n(p l (t)n) (1.3) 

i=l 

The maximum of absolute value of this entropy is equal to — Inn. So we 
normalize as follows: 

= (1-4) 

In section 0] we shall consider the hard disks systems. We shall compute 
an entropy functional similar to the space-entropy on extended torus with 
several cells. The probabilities are defined as for the space entropy in the 
Lorentz gas. We shall also do some comparisons of the if -theorem with the 
sum of normalized positive Lyapounov exponents. 



2 Entropy for collision map 

The entropy for the collision map is computed for a beam of N particles on 
a toric checkerboard with n cells. We start to calculate the entropy, just 
after all particles have executed the first collision. In this computation, all 
particles have the same initial velocity and are distributed in a small part 
of one cell. For each particle we determine the first obstacle and the angles 
(Pi, ipi) of the velocity incoming vector Vi(Pi) ( see the figures given in the 
appendix). For a uniform partition V of the space of the variables (P,ip), 
the entropy <S(t)is computed iteratively just after all particles have executed 
the t th collision. We use the formula where 

rPi+i ri>i+\ 

a(Vi) = / cos pdpdtp (2.5) 

is the invariant measure jTE] of the cell Vi = [Pi, A+i[ x bPi, ^i+ii an d Vt^Pi) 
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Figure 2: Entropy of the collision map versus number of collisions for (a) a beam of 
640 particles for a radius a=0.2, neighboring disks centers distance 1 and a partition of 
(f3,ip) space into 25 x 25 cells, (b) a beam of 512 particles for the obstacles of radius 0.2, 
neighboring disks centers distance 1 and a partition of (/?, tp) space into 9x9 cells. 



is the probability that a particle is located after t collisions in Vi computed 

as 

Number of particles in V\ having made t collisions 

N ' 
The velocity after the collision is computed from the following equation: 

V(P 2 ) = V(Pi) - 2(V(Pi).n)n (2.6) 

where n is the normal vector at the collision point. We explain in the ap- 
pendix the main geometric formula used for this computation. This entropy 
increase is shown in the Fig. |2]for various partitions and various initial distri- 
butions. The absolute value of the entropy of a distribution of particles, that 
we call its amount of entropy, represents in fact its distance to equilibrium. 
This is illustrated in the examples of randomly distributed initial velocity of 
particles having small amount of entropy (see Fig. |3] ) comparatively with 
beams of particles. It is to be noted that the amount of entropy increase 
under one collision is remarkably greater for the few first ones (more or less 
2-4 collisions) which corresponds to an exponential type increase (Fig. HJ. 
In order to calculate Lyapounov exponents by using the method of Benettin 
et al 0, first we calculate the Jacobian matrix in the tangent space of the 
collision map: 
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Number of collisions 
(a) 



Number of collisions 
(b) 



Figure 3: (a) and (b) are the entropy of the collision map with random initial conditions 
versus number of collisions for the system of particles of the Fig. |2 respectively. 




Figure 4: Logarithm of the collision map entropy versus number of collisions for the 
system of particles of the Fig. 



d/3i dcti 



\ dct2 da 2 J 

Now, comparing A S = max( l 5(t+ 1) — S(t)) (where the "max" is taken over 
t) with the positive Lyapounov exponent, A, of the collision map we verify 
the inequality: 

A S < X (2.7) 



as shown in Fig. |5J where this exponent is ~ 3.2. The maximal entropy 
increase by collision for the distribution computed in this figure is not far 
from this value. So it could be conjectured that in some suitable refinement 
limit, the entropy increase of a beam tends to the positive Lyapounov expo- 
nent. The rate of the approach to equilibrium is thus related to the positive 
Lyapounov exponent. Furthermore, the value of Lyapounov exponent is only 
dependent of — , i.e. the ratio of the distance between two successive obsta- 
cles over the radius of the obstacle, and its variation is exponential as shown 
in Fig. 

In order to compare the entropy increase as a function of the collisions 




Figure 5: (a) Lyapounov exponent and (b) entropy of the collision map, versus of number 
of collisions for each particle. We see that the maximum of the entropy increase between 
two collisions is less than of the value of the Lyapounov exponent. 



with the entropy increase as a function of time, we compute the distribution 
of mean free time for the first 3 collisions. From time histogram for the 
first three collisions of this system ( Fig. |7|), we see that a great number of 
particles have the same mean free time. As shown in the table ^ the mean 
free time vary during the first three or four collisions but after those, for the 
following collisions, rapidly the system comes near the equilibrium, where we 
have a constant mean free time approximately. 
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Figure 6: Lyapounov exponent versus D/a. 



Collision number 


1 


2 


3 


4 


5 


6 


7 


Mean free time 


1.966 


26.174 


5.801 


3.820 


3.611 


4.452 


4.177 


Collision number 


8 


9 


10 


11 


12 


13 


14 


Mean free time 


4.162 


4.208 


4.212 


3.863 


4.272 


4.051 


4.397 



Table 1: Mean free time obtained for a beam of 640 particles for a radius a=0.2, neigh- 
boring disks centers distance 1 and a partition of (/3, -0) space 25 x 25 cells. 



(a) 

Figure 7: Free time histogram for (a) first, (b) second and (c) third collision. 




3 Spatially extended Lorentz gas entropy 

The computation of the normalized space entropy equation by using (jl.4j) 
versus the number of collisions shows a remarkable exponential increase both 
for beams and for random initial distributions (Fig. |Hj). The computation of 
sum of the two positive Lyapounov exponents of the flow of one particle is 
equal to 1.046. Thus, we observe that the inequality between the normalized 
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increase of the density of the space entropy and this sum is verified. 




Figure 8: (a) Normalized space entropy of the Lorentz gas versus number of collisions for 
a beam of 640 particles for obstacles of radius a=0.2, neighboring disks centers distance 1 
and a partition of (x, y) space into 25 x 25 cells, (b) Logarithm of the space entropy versus 
number of collisions for this system. 



4 Hard disks 

Considering a uniform space partition of a large toric space we compute 
the particles densities, p iy and the normalized space entropy as a function of 
time by using the equation (jl.4|) . Starting with a distribution of disks with 
localized positions in some cell and random velocities, we compute binary 
collisions instants and the trajectories of the hard disks. These instants are 
determined by checking the distance between particles, after a time interval 
is passed. The Lyapounov exponents of the flow are calculated by using the 
Benettin et al. algorithm. The result is shown in the Figs. Eland E3 These 
figures show the entropy and logarithm of monotonic part of entropy versus 
time of the same gas with two distinct densities. The system in the Fig. 1101 is 
more dense than the system in FigEl and its entropy increases more rapidly. 
Fig. ^2 is a histogram of the number of collisions, so we see that the number 
of collisions in a fixed time interval is reduced for large time. From Figs. El 
and El we see that the monotonic part of the non-equilibrium entropy is also 
varying exponentially with respect to time. This shows that the collision is 
the main ingredient responsible of the entropy increase as described in the 
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Boltzmann equation theory. 

We shall now vary the density a = y and compute the characteris- 



Density 




As 


3.555 


0.367 


0.139 


0.889 


0.294 


0.115 


0.222 


0.239 


0.144 



Table 2: The data for the hard disks systems of radius, a = 0.05 and the same initial 
conditions, with cells 6 x 6, in terms of the density. 

tic quantities. The graph of the normalized positive Lyapounov exponents 
spectrum per particle for the same system as in Fig. El is shown in Fig 1 121 
The computation of the normalized sums of the positive Lyapounov expo- 
nent, jj J2\ t >o(~)T^~ )i shows that the inequality between maximum entropy 
increase and the sum of normalized of positive Lyapounov exponents is ver- 
ified ( Table El). 



s 




Figure 9: (a) Normalized space entropy and its monotonic part logarithm of the hard 
disks versus time for the 128 particles for the obstacles of radius a=0.05 which are initially 
localized in the first cell of (x, y) space with 6x6 cells and a density o\ — 0.889 disks per 
unit area. 
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(a) (b) 

Figure 10: Normalized space entropy and its and its monotonic part logarithm for the 
same system as Fig- with a density 02 = 3.555 disks per unit area. 
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Figure 11: Number of collisions histogram system versus time in Fig. El 

5 Concluding remarks 

The computations of the entropy amount of some given nonequilibrium 
initial distributions relatively to the equilibrium measure show an exponen- 
tial type increase for all considered partitions and distributions during initial 
stage after which the entropy increases slowly and fluctuates near its max- 
imal value. These computations confirm the existence of a relaxation time 
generally assumed in the derivation of kinetic equations 1J and the origin 
of the rapid increase of the entropy due to the number of collisions. The 
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Figure 12: Normalized spectrum Lyapounov of exponent of system in Fig. [5] 



dispersive nature of the obstacles is responsible of the exponential mixing 
type increase. This exponential type increase has been demonstrated for 
the Sinai entropy functional jTH] m hyperbolic automorphisms of the torus. 
On the other hand, the relation of the entropy increase to Lyapounov ex- 
ponents can be understood through Pesin relation and Ruelle inequality. In 
fact, the rate of entropy increase should be bounded by the Kolmogorov- 
Sinai entropy and such bound have been found by Goldstein and Penrose 
for measure-theoretical dynamical systems under some assumptions |14j . An 
open question is to characterize the measures reaching the upper bound. 

Any entropy functional is not a completely monotonic function of time 
for any dynamical system. In order to define a completely monotonic entropy 
functional for a dynamical system some conditions on the dynamics should 
be imposed. We can first suppose the map T on a phase space X to be a 
Bernoulli system or, slightly more generally, a X-system. This means that 
there is an invariant measure \i and some partition £o of X such that T£o be- 
comes finer than £ ( we denote it: T£ > £o)- Using the notation: T n £ = £ n , 
we obtain a family of increasingly refined partitions, in the sense of the above 
order of the partitions. Moreover, £ n tends, as n — > oo, to the finest partition 
of X into points, and £ n tends, as n —>■ — oo, to the most coarse partition, 
into one set of measure 1 and another set of measure zero. A physical proto- 
type of a Bernoulli and a iT-system is the above billiard fHJ QUI . A geometric 
prototype of a Bernoulli and a i^-system is uniformly hyperbolic system with 
Sinai invariant measure ^7] . A non-equilibrium entropy for a family of ini- 
tial measures, using conditional expectations relatively to the partitions 
was first obtained as an equivalence between the unitary group evolution 
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and a semi-group of contraction operators in the space of square integrable 
functions L 2 (fi) successively for the baker transformation ^3], for Bernoulli 
systems [7j and for K-systems [T3]. Its extension to the space of measures 
in K-systems has been realized in [Bj. In differentiable hyperbolic dynamical 
systems where the fibers of the £ ra partitions are pieces of contracting fibers, 
the construction of such entropy functional results from a generalized coarse- 
graining with respect to these contracting fibers, each fiber being a piece of 
manifold of zero measure. 



6 Appendix 
6.1 Collision Map 

We shall give the formula of the collision map. We consider a particle which 
undergoes the first collision with the disk of center Oi with velocity ~Vi(pi) 
and the second collision with the disk of center O2 with velocity Vi(p2)- Two 
cases are possible. First, we consider non-crossing of the centers line as in 
the Fig. 1121 In this figure the angle PiP 2 M is a 2 — /3 2 = — ( a i — Pi), where 
M is such that MP 2 is parallel to O1O2. We can write 

PiM = P1P2 cos(/9i — «i) = d — a cos ot\ — a cos a 2 - (6-8) 

and 

P 2 M = P1P2 sin(/3i — ai) = a sin a,\ — a sin 02, (6.9) 
if we eliminate a 2 between these equations we arrive at 

(3 2 = arcsinf- sin(/5i — ai) +sin/3i]. (6.10) 
a 

In crossing case which we present in Fig. [Hwe see that the angle P2P1M is 
equal to a 2 — /?2 = a i ~ Pi, and the length of P2M is changed to: 

P 2 M = P1P2 sin(/?i — ax) = a sin a.\ + a sin a 2) (6-H) 

then, we have 

P2 = arcsinf- sin(/?i — a x ) — sin (3i\. (6-12) 

(X 

To obtain fa in the first collision between particle and obstacle Fig. El , we 
take d = OPi, j3\ = and ai = $ in the collision map. 
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Figure 14: crossing Collision. 




Figure 15: Particle obstacle Collision. 



6.2 Algorithm description 

In this section we describe the algorithm which we used in our program for 
Lorentz gas. We first define in the main of our program the initial conditions 
for the particles and the obstacles positions. In the second step, we compute 
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with which obstacle, a particle will collide: we measure the angle between 
velocity of particle and the line between this particle and the center of ob- 
stacle, OPi in Fig. EH if this angle is less than or equal to the angle between 
this line, OPi, and the tangent line on the circle, i.e. PiN, in brief if *& < (p 
in Fig. we have a collision. Now, we use the collision map equation (j6.10|) 
or fjCj. to obtain the collision angles, /3 2 , and a 2 (see Fig. Ho"]) . In this step, 
we can also obtain the length of arrow of our induced collision map, i.e. P1P2 
(see Fig. H3j). easily as: 

pa = QPi ~ a r a2 (6-i3) 

cosv 

where 0.1 = $2 — Then, we can calculate the time of flight of particle 




between two collisions, respectively, as t — PiP 2 /w. This provides the tra- 
jectory of a particles. 

Let us turn the computation of space entropy. When a particle arrives at a 
wall of the big torus, before it does a collision with an obstacle (see on the 
Fig. |TJ) trajectories are pursued until it undergoes a collision on the torus. 
We have to compute the position of the obstacle that the particle will hit 
(see Fig. IT7|) and the angle a,\ in the collision map, and to determine which 
type of collision, i.e. crossing or non-crossing case, will occur. We first find 
the angle of collision 

M0 2 

P2 = $ + V 9 = arcsinf sin#], (6.14) 

then we arrive at ati and «i as 
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J a™ — Pi + arcsin[^(sin/3 2 — sin/3i)], a 2 = P 2 + (Pi — «i) 
\ = /3i - arcsinf^sin/^ + sin/?i)], = P2 — (Pi - «i) 

(6.15) 

where the superscript "c" corresponds to crossing case and "n" corresponds 
to non-crossing case see equations (jfi.lOj) and ()6.12|) . respectively. In the 
above equations, the d parameter is unknown, and will be recognized it in 
the end of this appendix. If we subtract the above equations we obtain 

a\ — a" = — arcsin[-(sin/3 2 + sin/?!)] — arcsin [— (sin p 2 — sin Pi)] (6.16) 

We can see the above equation yields a\ — a™ < 0. It means that in the 
same conditions the angle a\ in the non-crossing is greater than crossing 
case. Also, we can get the same conclusion for a 2 , i.e. a 2 < a, 2 . Now, we 
initiate the algorithm in the non-crossing case and we find a\ and a 2 . If 
•d < a 2 , thus, we had a correct supposition, otherwise, we must consider the 
crossing case, and we re-calculate these angles. In order to find in this case 
the parameter d = | O1O2 1 , we calculate it by approximation method. The 
equation that recognize d is: 

d = vdt cos(Pi — ai) + a(cosai + cosa^) (6-17) 

where dt is the time of free flight of particle between two collisions, see Figs. 
(P and El). In the above equation we have two unknown variables, ot\ and 
a 2 . We use the zeroth approximation as 

dm vdt (6.18) 
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where we used a <C vdt. Now, we calculate the angles, ol\ and a 2) as men- 
tioned in above of this appendix. Then, we re-calculate d with the first 
approximation, and we can repeat this procedure. However, the convergence 
is very rapid. 
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